home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus Leser 15 / Amiga Plus Leser CD 15.iso / Tools / Development / yacas_alg / yacas_morphos / share / yacas / linalg.rep / code.ys < prev   
Encoding:
Text File  |  2002-03-13  |  6.0 KB  |  346 lines

  1.  
  2.  
  3. /* Binomials */
  4. 10 # Bin(n_IsNonNegativeInteger,m_IsNonNegativeInteger)_(m <= n) <-- n! /(m! *(n-m)!);
  5. 20 # Bin(n_IsInteger,m_IsInteger) <-- 0;
  6.  
  7. /* Levi-civita symbol */
  8. Function("LeviCivita",{indices})
  9. [
  10.   Local(i,j,length,left,right,factor);
  11.   length:=Length(indices);
  12.   factor:=1;
  13.  
  14.   For (j:=length,j>1,j--)
  15.   [
  16.     For(i:=1,i<j,i++)
  17.     [
  18.       left:=indices[i];
  19.       right:=indices[i+1];
  20.  
  21.       If (Equals(left,right),
  22.       [ factor := 0 ; ],
  23.       [
  24.         If(Not(Apply("<",{left,right})),
  25.         [
  26. /*
  27.           Swap(indices,i,i+1);
  28. */
  29.           indices:=Insert(Delete(indices,i),i+1,left);
  30.           factor:= -factor;
  31.         ]);
  32.       ]);
  33.     ];
  34.   ];
  35.   factor;
  36. ];
  37.  
  38. Function("Permutations",{result,list})
  39. [
  40.   If (Length(list) = 0,
  41.   [
  42.     result;
  43.   ],
  44.   [
  45.     Local(head);
  46.     Local(newresult);
  47.     Local(i);
  48.     head:=list[1];
  49.     newresult:={};
  50.     ForEach(item,result)
  51.     [
  52.       For(i:=Length(item)+1,i>0,i--)
  53.       [
  54.         DestructiveInsert(newresult,1,Insert(item,i,head));
  55.       ];
  56.     ];
  57.     newresult:=DestructiveReverse(newresult);
  58.     Permutations(newresult,Tail(list));
  59.   ]);
  60. ];
  61.  
  62.  
  63. Function("Permutations",{list})
  64. [
  65.   Permutations({{}},list);
  66. ];
  67.  
  68. Function("InProduct",{aLeft,aRight})
  69. [
  70.   Local(length);
  71.   length:=Length(aLeft);
  72.   Check(length = Length(aRight),"InProduct: error, vectors not of the same dimension");
  73.  
  74.   Local(result);
  75.   result:=0;
  76.   Local(i);
  77.   For(i:=1,i<=length,i++)
  78.   [
  79.     result := result + aLeft[i] * aRight[i];
  80.   ];
  81.   result;
  82. ];
  83.  
  84.  
  85. Function("CrossProduct",{aLeft,aRight})
  86. [
  87.   Local(length);
  88.   length:=Length(aLeft);
  89.   Check(length = 3,"OutProduct: error, vectors not of dimension 3");
  90.   Check(length = Length(aRight),"OutProduct: error, vectors not of the same dimension");
  91.  
  92.   Local(perms);
  93.   perms := Permutations({1,2,3});
  94.  
  95.   Local(result);
  96.   result:=ZeroVector(3);
  97.  
  98.   Local(term);
  99.   ForEach(term,perms)
  100.   [
  101.     result[ term[1] ] := result[ term[1] ] +
  102.       LeviCivita(term) * aLeft[ term[2] ] * aRight[ term[3] ] ;
  103.   ];
  104.   result;
  105. ];
  106.  
  107.  
  108. Function("ZeroVector",{n})
  109. [
  110.     Local(i,result);
  111.     result:={};
  112.     For(i:=1,i<=n,i++)
  113.     [
  114.       DestructiveInsert(result,1,0);
  115.     ];
  116.     result;
  117. ];
  118.  
  119.  
  120. Function("BaseVector",{row,n})
  121. [
  122.     Local(i,result);
  123.     result:=ZeroVector(n);
  124.     result[row] := 1;
  125.     result;
  126. ];
  127.  
  128. RandomIntegerVector(_count,_coefmin,_coefmax) <--
  129.   Table(MathFloor(coefmin+Random()*(coefmax+1-coefmin)),i,1,count,1);
  130.  
  131.  
  132. Function("Identity",{n})
  133. [
  134.     Local(i,result);
  135.     result:={};
  136.     For(i:=1,i<=n,i++)
  137.     [
  138.       DestructiveAppend(result,BaseVector(i,n));
  139.     ];
  140.     result;
  141. ];
  142.  
  143.  
  144.  
  145.  
  146. Function("DiagonalMatrix",{list})
  147. [
  148.   Local(result,i,n);
  149.   n:=Length(list);
  150.   result:=Identity(n);
  151.   For(i:=1,i<=n,i++)
  152.   [
  153.     result[i][i] := list[i];
  154.   ];
  155.   result;
  156. ];
  157.  
  158. Function("Normalize",{vector})
  159. [
  160.   Local(norm);
  161.   norm:=0;
  162.   ForEach(item,vector)
  163.   [
  164.     norm:=norm+item*item;
  165.   ];
  166.   (1/(norm^(1/2)))*vector;
  167. ];
  168.  
  169. Function("ZeroMatrix",{n,m})
  170. [
  171.   Local(i,result);
  172.   result:={};
  173.   For(i:=1,i<=n,i++)
  174.     DestructiveInsert(result,i,ZeroVector(m));
  175.   result;
  176. ];
  177.  
  178.  
  179. Function("Transpose",{matrix})
  180. [
  181.   Local(i,j,result);
  182.   result:=ZeroMatrix(Length(matrix[1]),Length(matrix));
  183.   For(i:=1,i<=Length(matrix),i++)
  184.     For(j:=1,j<=Length(matrix[1]),j++)
  185.       result[j][i]:=matrix[i][j];
  186.   result;
  187. ];
  188.  
  189.  
  190.  
  191. /* Cramer-like matrix operations */
  192. Function("Determinant",{matrix})
  193. [
  194.   Local(perms,indices,result);
  195.   indices:=Table(i,i,1,Length(matrix),1);
  196.   perms:=Permutations(indices);
  197.   result:=0;
  198.   ForEach(item,perms)
  199.      result:=result+Factorize(i,1,Length(matrix),matrix[i][item[i] ])*
  200.                     LeviCivita(item);
  201.   result;
  202. ];
  203.  
  204. Function("CoFactor",{matrix,ii,jj})
  205. [
  206.   Local(perms,indices,result);
  207.   indices:=Table(i,i,1,Length(matrix),1);
  208.   perms:=Permutations(indices);
  209.   result:=0;
  210.   ForEach(item,perms)
  211.      If(item[ii] = jj,
  212.        result:=result+
  213.          Factorize(i,1,Length(matrix),
  214.          If(ii=i,1,matrix[i][item[i] ])
  215.                   )*LeviCivita(item));
  216.   result;
  217. ];
  218.  
  219.  
  220.  
  221. Minor(matrix,i,j) := CoFactor(matrix,i,j)/(-1^(i+j));
  222.  
  223. Function("Inverse",{matrix})
  224. [
  225.   Local(perms,indices,inv,det,n);
  226.   n:=Length(matrix);
  227.   indices:=Table(i,i,1,n,1);
  228.   perms:=Permutations(indices);
  229.   inv:=ZeroMatrix(n,n);
  230.   det:=0;
  231.   ForEach(item,perms)
  232.   [
  233.     Local(i,lc);
  234.     lc := LeviCivita(item);
  235.     det:=det+Factorize(i,1,n,matrix[i][item[i] ])* lc;
  236.     For(i:=1,i<=n,i++)
  237.         [
  238.          inv[item[i] ][i] := inv[item[i] ][i]+
  239.            Factorize(j,1,n,
  240.              If(j=i,1,matrix[j][item[j] ]))*lc;
  241.         ];
  242.   ];
  243.   Check(det != 0, "Zero determinant");
  244.   (1/det)*inv;
  245. ];
  246.  
  247. Function("Trace",{matrix})
  248. [
  249.   Local(i,result);
  250.   result:=0;
  251.   For(i:=1,i<=Length(matrix),i++)
  252.     result:=result+matrix[i][i];
  253.   result;
  254. ];
  255.  
  256. x X y := CrossProduct(x,y);
  257. x . y := InProduct(x,y);
  258. Commutator(x,y):=x*y-y*x;
  259.  
  260.  
  261. /* SylvesterMatrix */
  262.  
  263. Function("SylvesterMatrix",{poly1, poly2, var})
  264. [
  265.   Local(i,m,p,q,y,z,result);
  266.   y:=Degree(poly1,var);
  267.   z:=Degree(poly2,var);
  268.   m:=y+z;
  269.   p:={};
  270.   q:={};
  271.   result:=ZeroMatrix(m,m);
  272.  
  273.   For(i:=y,i>=0,i--)
  274.     DestructiveAppend(p,Coef(poly1,var,i));
  275.   For(i:=z,i>=0,i--)
  276.     DestructiveAppend(q,Coef(poly2,var,i));
  277.  
  278.   For(i:=1,i<=z,i++)
  279.   [
  280.     Local(j,k);
  281.         k:=1;
  282.     For(j:=i,k<=Length(p),j++)
  283.         [
  284.           result[i][j]:=p[k];
  285.           k++;
  286.         ];
  287.   ];
  288.  
  289.   For(i:=1,i<=y,i++)
  290.   [
  291.     Local(j,k);
  292.         k:=1;
  293.     For(j:=i,k<=Length(q),j++)
  294.         [
  295.           result[i+z][j]:=q[k];
  296.           k++;
  297.         ];
  298.   ];
  299.   result;
  300. ];
  301.  
  302.  
  303.  
  304. Function("MatrixRow",{matrix,row})
  305. [
  306.   Check(row > 0, "MatrixRow: row index out of range");
  307.   Check(row <= Length(matrix), "MatrixRow: row index out of range");
  308.  
  309.   Local(result);
  310.   result:=matrix[row];
  311.  
  312.   result;
  313. ];
  314.  
  315. Function("MatrixColumn",{matrix,col})
  316. [
  317.   Local(m);
  318.   m:=matrix[1];
  319.  
  320.   Check(col > 0, "MatrixColumn: column index out of range");
  321.   Check(col <= Length(m), "MatrixColumn: column index out of range");
  322.  
  323.   Local(i,result);
  324.   result:={};
  325.   For(i:=1,i<=Length(matrix),i++)
  326.     DestructiveAppend(result,matrix[i][col]);
  327.  
  328.   result;
  329. ];
  330.  
  331. Function("GenMatrix",{func,m,n})
  332. [
  333.   Local(i,j,result);
  334.   result:=ZeroMatrix(m,n);
  335.  
  336.   For(i:=1,i<=m,i++)
  337.     For(j:=1,j<=n,j++)
  338.           result[i][j]:=Apply(func,{i,j});
  339.   
  340.   result;
  341. ];
  342. HoldArg("GenMatrix",func);
  343.  
  344.  
  345.  
  346.